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I. INTRODUCTION 



We derive a tidal potential for a self-gravitating fluid star orbiting Kerr black hole along a timelike 
geodesic extending previous works by Fishbone and Marck. In this paper, the tidal potential is 
calculated up to the third and fourth-order terms in R/r, where R is the stellar radius and r the 
orbital separation, in the Fermi-normal coordinate system following the framework developed by 
Manasse and Misner. The new formulation is applied for determining the tidal disruption limit 
(Roche limit) of corotating Newtonian stars in circular orbits moving on the equatorial plane of 
' Kerr black holes. It is demonstrated that the third and fourth-order terms quantitatively play an 

important role in the Roche limit for close orbits with R/r > 0.1. It is also indicated that the Roche 
£S| ■ limit of neutron stars orbiting a stellar-mass black hole near the innermost stable circular orbit may 

^ ' depend sensitively on the equation of state of the neutron star. 

rjP '. 04.25.Dm, 04.25.-g, 04.40.Dg, 04.70.-s 
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Black hole is the most compact object in the universe and can tidally disrupt ordinary stars and compact stars 
> such as white dwarfs and neutron stars. A star of mass m and radius R will be tidally disrupted by a black hole of 
,— 1 ' mass M for the case 

o : 

m . m f r \ 3 

o 

O" 1 where r denotes the orbital separation and /i cr it is a constant sw 10-20 which depends on r/M, R/M, spin of the 
$_i , black hole, and equations of state (see Sec. V). If the condition (1) is satisfied outside a minimum orbital radius (e.g., 
60- r = QGM/c 2 for a star in a circular orbit around a Schwarzschild black hole where c and G are the speed of light and 
gravitational constant), a star will be tidally disrupted. According to Eq. (1), (i) an ordinary star of mass ~ M Q 
plunging inside a tidal radius of a supermassive black hole of mass smaller than ~ 10 8 M Q will be tidally disrupted, 
(ii) a white dwarf of mass ~ 0.7Mq and radius ~ 10 4 km plunging into a massive black hole of mass smaller than 
~ 10 5 Mq will be tidally disrupted, and (iii) an inspiraling neutron star of mass ~ XAMq and radius ~ 10 km in 
a circular orbit around a stellar-mass black hole of mass smaller than ~ 5M will be tidally disrupted before the 
neutron star reaches the innermost stable circular orbit (ISCO). During the tidal interaction of the ordinary star by 
a supermassive black hole, an ultraviolet flare of a characteristic light-curve may be emitted at the center of galaxy 
and be observed (e.g., [1,2] for a review). White dwarfs or neutron stars tidally disrupted by a stellar-mass black hole 
will form a massive disk which may be a possible candidate of the central engine of gamma-ray bursts [3] . Detection 
of gravitational waves at a tidal disruption of neutron stars by a stellar-mass black hole may constrain the equation 
of state of neutron stars [4,5]. These examples show the tidal disruption of ordinary stars and compact objects by a 
black hole is likely to happen frequently and be an interesting phenomenon in the universe. This fact stimulates the 
theoretical study for the tidal disruption of stars by a black hole. In this paper, we present a new general relativistic 
formulation with higher-order corrections of the tidal potential which can be used to clarify the criterion of the tidal 
disruption for close orbits more accurately than that by previous works. 

Numerical studies for the tidal disruption of a star by a black hole have been extensively performed assuming 
the Newtonian (e.g., [6-11]) or post-Newtonian [1] gravity with a point particle approximation for the black hole. 
However, such an approximation is not quantitatively appropriate for analyzing the tidal disruption near the black 
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hole, since general relativistic effects are essential for such close orbits. In [12], a numerical result for a general 
relativistic simulation was presented, but the authors ignored the self-gravity of the fluid star assuming that it is 
much weaker than the tidal force of the black hole. Relativistic tidal problems have been widely studied in the 
so-called tidal approximation (e.g., [13-17]). In this approximation, one assumes that the mass of a star m is much 
smaller than the black hole mass M and that the stellar radius R is smaller than the orbital radius r (or the curvature 
radius of the black hole spacetime). As a result of these assumptions, one may assume that (i) the center of mass of 
the star moves around the black hole along a timelike geodesic in the black hole spacetime and (ii) the tidal field from 
the black hole is computed from the Riemann tensor of the black hole spacetime in terms of the geodesic deviation 
equation. In addition to (i) and (ii), one often assumes that the self- gravity of the star is described by the Newtonian 
gravity. This approach has been used for studies of the tidal disruption limit of ordinary stars and compact objects 
[9,11,13,17], and for hydrodynamic simulations of the tidal disruption of ordinary stars and white dwarfs [16]. The 
purpose of this paper is to improve a calculation of the tidal potential in this framework, since the improvement is 
necessary for some problems as explained in the following. 

In the tidal approximation, one assumes R <C r and then expands the tidal potential in terms of R/r. This 
approximation is illustrated for the Newtonian potential <j> as follows: Considering of a point particle of mass M p 
at (r, 0, 0) and expanding it around the origin, one obtains 

GM p 



\J{x - r) 2 +y 2 + z 2 



GM p 



x 2x 2 -y 2 - z 2 x(2x 2 - 3y 2 - 3z 2 
1-| 1 2 1- _i % 

r 2r 2 2r 3 



8x 4 - 24x 2 (y 2 + z 2 ) + 3y 4 + 3z 4 + 6y 2 z 2 _ 5 . 

8r 4 + U ^ r ) 



(2) 



where we assume that \x\, \y\, \z\ <C r (i.e., R <C r). In the standard tidal approximation, one takes into account terms 
up to the second order in R/r neglecting the terms of 0[(R/r) 3 ], resulting in 



^tidal approximation 



GM p 



■A ~ - I i ~ 

r 2r 2 



x 2x 2 — y 



(3) 



The standard tidal approximation works for determining the tidal disruption limit in many problems, but docs not 
in some problems: Equation (1) implies that a tidal disruption occurs if the following condition is satisfied; 

7 £ V". (4) 

Here, Q = m/M . Thus, with increasing the mass ratio Q, the critical value of R/r for the tidal disruption increases: 
For Q > 1CP 3 , the tidal disruption sets in for R/r > 0.05, indicating that neglecting more than third-order terms in 
R/r yields an error of > 5% and may not be a very good approximation for the computation of the tidal disruption 
limit. For example, the tidal disruption of white dwarfs by an intermediate mass black hole of ~ 10 3 M Q is the 
case. Even for Q < 10~ 3 , the star will be elongated during the close-encounter with a black hole in parabolic or 
highly elliptical orbits. In such cases, R will increase with decreasing r, and hence, the higher-order terms may be 
important. For binaries of a black hole of mass smaller than ~ 5M Q and a neutron star in close quasicircular orbits, 
the higher-order terms are also important since R/r > 0.2 at the ISCO of radius ~ 6GM/c 2 . Our numerical results in 
the framework of Newtonian gravity indeed suggest that the third- and fourth-order terms play an important role for 
a system of R/r > 0.1 [18]. For a rapidly rotating black hole, the radius of the ISCO can be as small as ~ GM/c 2 . In 
this case, the higher-order terms in R/r may become quite important. These facts illustrate that it is often necessary 
to take into account the higher-order corrections of the tidal potential for computation of quantitatively improved 
results in the tidal problem. 

In this paper, we derive a general relativistic tidal potential induced by a black hole in which higher-order corrections 
are taken into account. In contrast to the previous works [13,15], we do not use the geodesic deviation equation to 
calculate the tidal potential. Instead, we derive the tidal metric for an observer moving along a timelike geodesic on 
a black hole spacetime in the Fermi normal coordinates following Manasse and Misner [19]. In this method, the tidal 
potential can be calculated from the tidal metric in a straightforward manner. 

Using the new formulation, we numerically compute equilibrium states and a tidal disruption limit (Roche limit) 1 
for corotating stars of polytropic equations of state in a close circular orbit around a Kerr black hole. Comparing the 



The Roche limit is defined as the tidal disruption limit for a star of corotating velocity fields. 



2 



numerical results with those computed by the standard tidal approximation, we illustrate that the higher-order terms 
of the tidal potential play a quantitatively important role for R/r = 0(0.1). The numerical results in this framework 
will be also used to calibrate accuracy of a numerical result of fully general rclativistic quasicquilibrium states for a 
binary of a black hole and a neutron star which will be computed in near future. (It should be noted that there has 
been no comprehensive work about the tidal disruption of stars by a black hole in full general relativity, although 
there are primitive works for binaries of a black hole and a neutron star [21,22].) The solution obtained in this work is 
valid for r ^> R and M 3> m, and hence, the calibration can be carried out for the case of low-mass neutron stars. In 
addition, dependence of the tidal disruption limit on the equations of state for the star is investigated. It is indicated 
that the tidal disruption limit of a neutron star by a stellar-mass black hole depends sensitively on the equation of 
state of the neutron star. 

The paper is organized as follows. In Sec. II, we derive a general expression of the metric for an observer moving 
along timelike geodesies in the Fermi normal coordinates. In Sec. Ill, the results of Sec. II are applied to the Kerr 
spacetimc. In Sec. IV, the tidal potentials are computed for equatorial circular orbits. In Sec. V, the tidal disruption 
limit (Roche limit) of corotating stars of equatorial circular orbits around a Kerr black hole is presented for a wide 
range of the spin parameter of the black hole and the equations of state for the star. Sec. VI is devoted to a summary. 

In the following, we adopt the geometrical units c = G = 1. Latin and Greek indices denote spatial and spacetime 
components with r — x° as the time coordinate in the Fermi normal frame, g^, T^ a , and R^^s denote the spacetime 
metric, Christoffel symbols, and Riemann tensor, respectively. 6ij(= 5^) denotes the Kronecker delta. The comma 
and semi-colon denote the ordinary and covariant derivatives. For simplicity, we often use the notations [23] 

A(ijk) — g(^ijfc + Ajik + Ajki + Akji + Akij + Aikj), (6) 

A( ai ... ai ) = 77 y]Ai/(i)-a/ f i)i (7) 

L i 

where in the last equation, the sum is taken for all permutations. 

II. EXTERNAL METRIC IN THE FERMI NORMAL COORDINATES 

A. Outlines and definition of the Fermi normal coordinates 

We derive the metric in an observer frame which moves along timelike geodesies around a Kerr black hole. As 
Manasse and Misner [19] show, this can be done by finding the relation between the Riemann tensor and the metric 
in the Fermi normal coordinate system which is one of the local inertial frames. Then, our goal is to write the metric 
in the neighborhood of an observer in this coordinate system as 

g^iiy ^StAViijK ^ ~\~ ~^Qiiv,ijk% 3^ ~t~ ~^^9nv,ijkl% ^ % % "i" 0{x ), (8) 

where rj^ is the flat metric and x l denotes a spatial coordinate in the Fermi normal coordinate system. The coefficients 
9tiv.ij--- are related to the Riemann tensor of the spacetime. Using the definition of the Riemann tensor and the geodesic 
deviation equations, Manasse and Misner [19] derived the quadratic term in Eq. (8). Our purpose is to derive the 
third and fourth terms following the method developed by them. 

Since the Fermi normal coordinate system is a local inertial frame, g^ v is the flat metric along the observer frame, 
and the Christoffel symbols vanish 

= o, (9) 
IV. = g^TZa = o. (10) 

In addition, the time direction in the Fermi normal coordinates is chosen as the direction of a timelike geodesic Q. 
Since Eqs. (9) and (10) are preserved along the timelike geodesic Q, one obtains 

Ffiva,0 = Ffiva,00 = Ffiva,0---0 = 0. (12) 
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From the definition of the Christoffel symbols, 

9fj,v,a — T/^^a "i" ^i/p,a: (13) 

one finds along Q, 

g^u.aO — g^iy,a00 = gpu,aO---0 = 0. (14) 

B. Relation obtained from the definition of the Rieraann tensor 

In this subsection, we present relations necessary for deriving the tidal potential up to the fourth order using the 
definition of the Riemann tensor 

r — r — r 4- v a r — r a r (^ K\ 

L ^pvpa — 1 crpp.v L ovp.p i x vp L apa 1 ^p 1 - Qi/a- \ ±0 J 

In the following, the components of the Riemann tensor along Q are computed. For fi = in Eq. (15), we find that 
the nonzero components are 

Roiuj = = — rVojji; (16) 

where Eqs. (9)-(12) are used. For v = in Eq. (16), one finds 

500,ij = — 2r :) 00,i = — 2-RoiOj": (17) 

where Eq. (14) is used. This gives one of the second-order terms in Eq. (8). 

Since the Christoffel symbols are vanishing, the covariant derivative of the Riemann tensor is equal to the ordinary 
derivative along Q as 

Rpvpa\a Rpupa.a- (1$) 

Thus, from Eq. (15) and T^ p = 0, one obtains a relation along Q 

Rp.vp<7;a r \jpp,vcx ^ aup.pa: (19) 

and thus, 



Ri0j0;k — Foijfik — FoOj,ik 
1 

2 



— o ( 9oi,jko + <7cy,ifco - goo,ijk ) • (20) 



Equation (20) leads to 



3 

RoiOj-k + RojOk:i + RokOi:j = {90i,jk + goj.ki + g0k,ij),0 — ^900,ijk- (21) 



In Sec. II D, this equation is used for deriving the third-order terms in Eq. (8). 
The second covariant derivative of the Riemann tensor along Q is written as 

Thus, 

Rpvpa:(af3) = Rpupa.af) — r^^ a j Ryi/pa ~ ^Z((3,a) Rp-lpa ~ ^p(/3,a)-^M I/ 7<^ — ^ct(/3,q) Rp-vpl ' (23) 

and 

RoiOj;(ki) = RoiOjM - ^1(k,i) R 'rOjo - r^ fc ^(Ri-yjo + Rj~/io) - FJ^R-yOio- (24) 
Using Eq. (15), one obtains 
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and hence, 



RoiOj,kl — ^ (^90i,0jkl + gOjfiikl — gOO,ijkl ~ 9ij,kW0 

+^Oi,k^aOj,l + ^M,l^aOj,k — ^O0,l^o:ij,k — ^QO,k^ocij,U (25) 

RoiOj,kl + RoiOkJl + RoiOLjk + R()jOk,il + RojOl,ik + RokOLij 

3 

— 2 y90i,jkl + SOj.ifci + <70fc,ij'i + 90l,ijk),0 - ^900,ijkl 

1, , 

- ^ + 5ife,ji + ftijfc + ffjfc.ii + ffji.ifc + 9kl,ij),00 

+ 4 ( r 0( l j) r a0(fcJ) + r 0(i,fc)r Q 0(j,/) + r 0(iJ) rQ °(j'^)) 

- 3 ( r oo,i r a(j/s,o + roo,jr a (j feii ) + fc r Q ( lii() + r^r^^)). (26) 

In Sec. II D, this equation is used to derive the fourth-order terms in Eq. (8). 

C. Relations obtained from the geodesic deviation equations 

In this subsection, relations necessary for calculating the tidal potential up to the fourth order are derived using 
the geodesic deviation equations. 

For the tangent of a geodesic, u^, and the displacement vector to an infinitesimally nearby geodesic the geodesic 
deviation equation is written as 

/V,(/V„/) = -R^uWz". (27) 

Note here that the geodesic that we consider is not restricted to Q. Using an affine parameter A for the geodesic (i.e., 
u» = {d/dXY), Eq. (27) is rewritten to 

+ 2^ro CT + (r<^ CT + r^r^ - r^r^)«"«^ Q = -i? raa <w A (28) 

or 

+ 2^r(>- + T^ a u^z a = 0. (29) 

Now, we consider the family of spacelike geodesies in the Fermi normal coordinates of the form 

x° = const, and x 1 = a l \ (30) 
where a 1 is a constant three-component. Using its definition, is written to 

(31) 

As the displacement vector z° ', we choose the spatial vector defined by 

.s?) = A5 "- (32) 

Substituting Eqs. (31) and (32) into Eq. (29) leads to 

2Y$ j ot + XT%^a k = 0. (33) 

Our purpose here is to derive the relations among the derivatives of the Christoffel symbols along Q which will be 
useful in the subsequent calculations. To obtain the relations, Eq. (33) should be evaluated for A = 0. However, Eq. 
(33) is trivial because of the vanishing Christoffel symbols in the Fermi normal coordinates. Thus, we carry out a 
Taylor expansion of Eq. (33) assuming that A is small, and the first-, second-, and third-order terms in A provide 
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2rr jlfc ^ fe +r^.a^ fe = 0, (34) 

rf^aW+T^aW^O, (35) 

|r&, fcIn aWa» + ^ Mi „aWa" = 0. (36) 

Equations (34)-(36) lead to cyclic relations among the derivatives of the Christoffel symbols along Q as 

r^y^) = 0, (37) 

^n(ij,kl) = °> ( 38 ) 

^fi{ij,kln) = 0- (39) 
Setting ^ = in Eqs. (37)-(39), the following useful relations are derived; 

9o(i,jk) = 0, (40) 

3{goi,jkl + 90j,ikl + gokdjl + 9M,ijk) = {9ij,kl + 9ik,jl + 9il,jk + 9jkM + 9jl,ik + 9kl,ij),0, (41) 
HgOijkln + 90jAkln + 90kajln + 90l,i]kn + 90nAjkl) 

{9ij-kln ~t~ 9ik,jln ~t~ 9il,jkn 9in,jkl "i" 9jk,iln "i" 9jl-ikri "i" 9jn.ikl ~t~ 9kl,ijn ~t~ 9kn,ijl ~t~ 9ln,ijk) ,0 • (42) 

Using Eq. (28) with the same strategy as that for deriving Eqs. (34)-(36), one can also obtain relations among the 
Christoffel symbols and Ricmann tensor as 

V^^r^a*, (43) 

^•fcV afca, = 2r £-,fci aiafca '' ( 44 ) 

V"*^" = fr£. w „<* J Wa" + (r?. ,1^ - 1^1^ JaWa". (45) 

Here, rj fe n a : 'a' c Q:™ = because of Eq. (37). Thus, the last term of Eq. (45) vanishes. 

From Eqs. (43)-(45) together with Eqs. (37)-(39), one obtains the relations between the Christoffel symbols and 
the Riemann tensor along Q as 



■ ij,k^ x ik,j ~ 3 I Rjik + ^kij ^ ) ' (^6) 



^ijM + rffcjz + rfi jk — 4 (j^jik M ,i + ^jii ^fe + ^feij + ^/ij ^fe + ^Zifc *J + ^kil ^jj ' 



(47) 



ij,kln "r" ik,jln * il,jkn * injkl 



r I i(?fe) .In .fen i in ,fc( i(fcf) .in i(fen) ,i( 



2 

5 I ± ^'^{jk) -In r ±tj i(jl) ,kn ^ ±l, i(jn) ,kl T ±>J i(kl) ,jn T JL i(fen) ,j7 ^ ±l, i(ln) Jk 

I r u 1 r" r'' 1 r" r' : 1 r" 1-'' 1 r" r'' r ,; ' r'' ' ! < v. 

~~ 5 I L i(j,k) L u(l,n) + 1 i(j,l) L u(k,n) + L i(j,n) L u(k,l) + 1 i(k,l) L v (j,n) + L i(k,n) L + 1 i(l,n) L v {j,k) I M:1 - 1 



Using Eq. (37), Eq. (46) is rewritten to 



3 

For fi = in Eq. (49), 



2 

r w -fe,i = -^Ri(jk)n- (49) 



2 1 1 

-^Ri(jk)a = ^Ojk,i = - (goj.ik + 90k,ij) = —-^9$i,jk-, (50) 



where Eqs. (40) and io = are used. Thus, 



3 

This gives one of the second-order terms in Eq. (8). 



gok,ij = n ( Roijk + Rojik ) • (51) 
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Equation (49) is also used to derive the relation for gui,ij along Q 



9kL 



Tfciij + T/fcij — g ( Riklj + Rilkj ) • 



(52) 



This also gives one of the second-order terms in Eq. (8). 
From Eqs. (52), one finds the following relations along Q; 

9i(j,kl) = 0, 
9(kl,j)i = 0. 

Equations (40), (53), and (54) are preserved along Q, i.e., 

\3o(k,ij)]fi = 0, 
\9i(j,kl)],0 = 0, 
[9(kl,j)i},0 = 0. 

Substituting the last two relations into Eq. (41), one obtains 

90(i,jkl) = 0. 

These relations are useful for deriving the third-order terms in Eq. (8). 



(53) 
(54) 



(55) 
(56) 
(57) 



(58) 



D. Deriving the third- and fourth-order terms 

From Eqs. (21), (40), (47), (55), and (58), third derivatives of 00 and Oz components of the metric are found to be 

(59) 



g00,ijk — — g (^RoiOj:k + RojOk:i + RokOi;j^j — ~ 2-Ro(i|0| i 



goijkl — T ( RijkO;l + RikjO;l + RilkO;j + RiklO;j + RijlO;k + RiljO;k) — 7jRi(jk\0\;l) • (60) 



For ij components, 



9ij,kln — ^ijk,ln ~t~ ^-*jik,lr, 
1 



g [ ^ijkjn ~\~ r jil^kn ^ijn,kl ~t~ ^jik.ln ~\~ -f jil,kn ^jin : kl 



Rkjli-.n ~t~ Rkjni\l "i" Rljki;n ~t~ Rnjki;l ~t~ Rnjli;k "i" Rljni;k 



where we use Eqs. (38) and (47). 

Using Eqs. (24), (26), and (41), the fourth derivative of 00 component of the metric is written as 



goo,ijki 



1 



+ 3 



RoiOj;(kl) + RoiOk;(jl) + RoiOl;(jk) + RojOk;(il) + RojOl;(ik) + RokOl;(ij) 
RoiOj:(kl) + RoiOk;(jl) + RoiOh(jk) + RojOk;(il) + RojOl;(ik) + RokOl;(ij) 



R^iai\c,R 



R^lAl^nR,. 



- -1R0(i\Q\j-M) + 8R %i\0\ R \n\if)0- 

From Eq. (61) and R(ijk)i = 0, it is found 



(61) 



(62) 
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and hence, 



Using Eqs. (64) and (42), one obtains 



9(ij,kln) = 9i(j,kln) = 0, 
[9(ij,kln)],0 — [9i(j,kln)]fi = 0. 
90(i,jkln) = 0. 



By a straightforward calculation, one finds 

^fiij.kln ~\~ ^fj,ik,jln ~t~ r^/j/cn "i" ^fiinjkl — ^[9fJ-iJkln 9^(j,kln)i 9i(j,kln)fi} • 

For n = in Eq. (66), 

where Eqs. (63) and (65) are used. Thus, 



(63) 
(64) 
(65) 
(66) 

(67) 



9Qi,jkln — ^ ( *-0ij,kln ~t~ ^Oik ,jln ~t~ ,jkn ~\~ ^Qinjkl 



Y5 ( Ri(jk)0,ln + Ri{jl)0,kn + Ri(jn)0,kl + Ri(kl)0,jn + Ri(kn)0,jl + Ri(ln)0,jk 



4 

15 



_ (^i(j,k)^0u(l,n) + ^i(j,l)^0u(k,n) + ^i(j,n)^0v(k,l) 
+ ^i(k,l)^0u(j,n) + ^i(k,n)^0u(j,l) + ^i(l,n)^0v{j,k) 
Ri(jk)0,ln + Ri(jl)0.kn + Ri(jn)0M + Ri(kl)0,jn + Ri(kn)0jl + Ri(ln)0,jk 

Ri(jk) R®(ln)0 + Ri(jl) °Ro(kn)0 + Ri(jn) Ro(kl)0 + Ri(kl) Ro(jn)0 + Ri(kn) R®(jl)0 + Ri(ln) Ro(j k)0 

^ I D "it? I D "it? _1_P m R 

_ 135 \ i ^ fc > ^C") + K *0'0 H m(kn)0 + K i(jn) H m(kl)0 

_L P m P _L P m P _L P m P 1 



45 



- 8 B 16 R Op 16 fi> ™B 

- g-KiOfe|0|,;n) - Y5 K ^ k K Win)0 ~ ^^i(jfc R\m\ln)0, 

where we use Eq. (48). 

The fourth derivative of ij components is derived as follows: 

9ij.klmn — ^ijk,lmn ~t~ ^jik,lmn 
1 



(68) 



Ri(kl)j.mn ~t~ Ri(km)j,ln ~t~ Ri(kn)j ,1m ~t~ Ri(lm)j.kn Riiln)j,km ~t~ Ri(mn)j,kl 



^t(k,l)^ jv(m,n) +^i(k,m)^ji'(l,n) +^i(k,n)^j^(l,m) 
+ r i(i,m) r ^(fc,«) + r i(i,n) r w(fe,m) + r j(m,n) r ^(fc,0 



Ri(kl)j,mn ~t~ Ri(km)j,ln ~t~ Ri(kn)j.ln Ri(lm)j.kn ~t~ Ri(ln)j,km Ri(mn)j. 



kl 





00 components 


equation 


Oi components 


equation 


ij components 


equation 


2nd order 


500, ij 


(17) 


90i,jk 


(51) 


gij,ki 


(52) 


3rd order 


gOO.ijk 


(59) 


gOiJkl 


(60) 


gij.klm 


(61) 


4th order 


S00,ijkl 


(62) 


g0i,jklm 


(68) 


gij,klrnn 


(69) 



TABLE I. Equation numbers from which one finds the relation between g^v^jk— an d the corresponding Riemann tensor in 



the Fermi normal coordinates. 



2 

15 



2 

"45 



R i(kl)° R j(mn)0 + R i(km) R j(ln)0 + R i(kn) R j(lm)0 
+ R i(lm) R j(kn)0 + Ri(l n ) R j(km)0 + R i( mn ) R j(kl)0 
R i(kl) ,R j(mn)p + R i(km) R 0(ln)p + R i(kn) R i(l™)P 



+ R i(lm) R j(kn)p + R i(ln) PR j(km)p + R i(mn) R i(kl)p 



~R< 



(kl\j\, 



^iikl - n -b'|mn)0 "|j|mn)p- 



(69) 



We note that the second ordinary derivatives of the Riemann tensor in Eqs. (68) and (69) are transformed to the 
covariant derivatives using 



i(jk\0\;ln) ~ R i(jk\0\,ln) ~ ^ n i(nl' n \<r\jk)0, 
R i(jk\m\;ln) R i(jk\m\,ln) ~n R i{nl R \u\jk)m- 



(70) 
(71) 



This implies that goijkim and gij t klmn can be written in the covariant form. 

To summarize this section, we derive the coefficients g^ v ^j... in Eq. (8) which denotes the metric at t = x° ^constant 
in the neighborhood of the timelike geodesic Q. This is used as the tidal field from a black hole for a star moving 
along Q. The equation numbers for the relations between the derivatives of the metric and the Riemann tensor in the 
Fermi normal coordinates arc summarized in Table I. 



III. COMPONENTS OF THE RIEMANN TENSOR FOR A KERR SPACETIME 

To compute components of the Riemann tensor and its covariant derivative for the Kerr metric in the Fermi normal 
coordinates, we adopt the method developed by Marck [15]. Namely, we first calculate the components in a standard 
tetrad frame of the Boyer-Lindquist coordinate system [24], and then, perform a coordinate transformation to the 
Fermi normal coordinates. 

The Kerr metric in the Boyer-Lindquist coordinate system is written as 

j2 / 2Mr\, , 4Mrasin 2 <9 (r 2 + a 2 ) 2 - Aa 2 sin 2 6 . 2 „ , 2 £, 2 ^ , n2 
ds 2 = - 1 — dt 2 dtdip + ± '—- sin 2 9d(p 2 + —dr 2 + Y,d9 2 , (72) 



A 

where 

S = r 2 + a 2 cos 2 6, A = r 2 + a 2 - 2Mr, (73) 
and M and a denote the mass and spin parameter. A standard tetrad for the Kerr metric is defined by 

(e (0) ) M = 0,0, -a sin 2 (74) 

(e (1) )^=(o,yf,0,o), (75) 
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(e (2) V= (o,0,VE,o) 

,(3h = f_ flSin " 



,0,0. 



(r 2 + a 2 ) sin 6* 



(76) 
(77) 



Note that the sign convention for (e^)^ is different from that of [15]. 

In the first step, we compute the components of the Riemann tensor and its covariant derivative in the standard 
tetrad. The nonvanishing components of the Riemann tensor in this tetrad frame are 



#(i)(2)(i)(2) - fl(i)(3)(i)(3) - 2 R{im{1){0) = ^2 R(2){3){2){3) 

Mr(3a 2 cos 2 9 - r 
E3 



-R 



(2)(0)(2)(0) 



-R 



(3)(0)(3)(0) 



1 aikfcos0(3r 2 - a 2 cos 2 9) 

«(1)(2)(3)(0) - -«(1)(3)(2)(0) - -2«(1)(0)(2)(3) = ^3 ■ 

The tetrad components associated with the first and second covariant derivatives are defined as 

Q(a)(b)(c)(d)(e) = ^ vRvpa\{e(a)T {e(b)Y {Z(C)Y (e {d) Y {e {e )) X , 

P(a)(b)(c)(d)(e)(f) = V( Q V ( 3 ) i? 1/(9aA (e( Q )) Q (e( b )) /3 (e (c) ) ,y (e (d) ) p (e (e) ) <J (e (/) ) A . 
The nonvanishing components of Q( a )(b)(c)(d)(e) are 

Q(o)(l)(2)(l)(2) = Q(o)(l)(3)(l)(3) = X<9(o)(l)(0)(l)(0) = -^Q(o)(2)(3)(2)(3) = -<9(o)(2)(0)(2)(0) = -Q(o)(3)(0)(3)(0) 



(78) 
(79) 



(80) 
(81) 



_ / 3MJ1A 1 / 2 121 



Ma rJi sin 9 cos 9 



-,o,o , 



Q(o)(l)(2)(l)(0) ~ Q(o)(2)(3)(3)(0) - 0,0, 



^ S9/2 ' £9/2 

YIMar J2A 1 / 2 cos 9 \2Ma 2 rJ 2 sin (9 cos 9 



£9/2 ' £9/2 

3MJ1A 1 / 2 3MaJisin<r 



/ 3MJ1A 1 / 2 3MaJisin6i\ 

<y(a)(l)(2)(2)(3) - -«(o)(l)(0)(3)(0) - ^ 0, ^ , ^ J 



Q(a)(l)(2)(3)(0) = -Q(o)(l)(3)(2)(0) = - Q(a)(l)(0)(2)(3) 



S 9/2 

12MarA 1 / 2 J 2 cos6» 3MaJisin(9 



£9/2 



£9/2 



,0,0,, 



, 3MaJ 1 sm9 YIMar J 2 A 1 / 2 cos 9 

Q(o)(l)(3)(l)(0) - -Q(a)(2)(3)(2)(0) - I ^ , ^773 > °, 



Q(o)(l)(3)(2)(3) ' Q(o)(l)(0)(2)(0) 



£9/2 

12Ma 2 rJ 2 sin0 cos 3MJ1A 1 / 2 



£9/2 



£9/2 



0,0 



where 



Ji=r 4 - 6a 2 r 2 cos 2 9 + a 4 cos 4 I 
J 2 = r 2 — a 2 cos 2 6*. 



(82) 
(83) 
(84) 
(85) 
(86) 
(87) 



(88) 
(89) 



For P( a )(b)(c)(d)(e)(f), the explicit form of the components in a general orbit is very complicated. Thus, we only write 
the nonvanishing components in the equatorial plane setting 9 = it/ 2; 



(o)(6)(l)(2)(l)(2) 



3M 

(o)(6)(3)(0)(3)(0) = — 



/ _(4 r 2 _ Q r M + 5a 2 ) 







3M 



(o)(6)(l)(2)(l)(3) = ^(o)(6)(2)(0)(3)(0) 



V 

f \ 

* A aA 1 / 2 

* * 

y * * * J 





A + 4a 2 

* 3A SaA 1 / 2 

* * -(Mr -3a 2 ) J 



(90) 



(91) 
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p p 3M 

^(o)(6)(l)(2)(l)(0) ~ ^(o)(6)(2)(3)(3)(0) - T^f 



/O 

* SaA 1 / 2 -(Mr -8a 2 ) 

* * 

* * * 



^(o)(6)(l)(2)(2)(3) ~ -^(o)(6)(l)(0)(3)(0) 



V 

3M 



/ 8A - Mr 8C1A 1 / 2 \ 



* 



V 

/ aA 1 / 2 a 2 \ 



P 



3M 



(o)(6)(l)(2)(2)(0) 



(o)(6)(l)(3)(3)(0) 



* 

* * 






J 



Y5 P (a)(6)(l)(2)(3)(0) - ~2j;-P(o)(6)(l)(3)(2)(0) 



1 



:P(. 



Pi 



(a)(6)(l)(3)(l)(3) 



-Pi 



(o)(6)(2)(0)(2)(0) 



(o)(6)(l)(3)(l)(0) 



= -Pi 



(o)(6)(l)(3)(2)(3) 



(o)(6)(l)(0)(l)(0) 



(a) (6) (2) (3) (2) (0) 



(o)(6)(l)(0)(2)(0) 



(a)(b)(2)(3)(2)(3) 



3M 



3M 

2T 7 



36 - (o)(6)(l)(0)(2)(3) ~ 

/ -(4r 2 -9rM + 7a 2 ) 



V 

/ IQaA 1 ' 2 







/ aAV2 

M * 

r 7 * * 

^ * * * J 



(3A + 4a 2 ) 

* A aA 1 / 2 

* * -Mr + a 2 I 

\ 



3M 
2^ 



6M 



* -16&A 1 / 2 

* * —Mr 
\ * * * J 

8r 2 - YJMr + 16a 2 

* 

* * 

* * * 

( -(4r 2 - 9Mr + 6a 2 ) 

* 2(A + 2a 2 ) 

* * 2A 
y * * * 



o 



2aA!/2 
-(Mr -2a 2 ) j 



(92) 



(93) 



(94) 



(95) 



(96) 



(97) 



(98) 



(99) 



Here, the components in the 2-matrix form of the subscripts (a) and (6) are shown in the order (1), (2), (3), and (0). 

To compute R a b c d, Q 'abode, and P a bcdef m the Fermi normal coordinates, we need to prepare the transformation 
matrix from the standard tetrad frame to the Fermi normal coordinate frame. Denoting it as A a a \ we have 

Rabcd = - R (a)(6)(c)(<i)Ai a) A['' ) A( c) A^ ) , 



labcde - Q(a)(h)(c)(d)(e)Ai a) A[ fc) A( c) A^ ) A^ e) , 



Pa 



bcdef 



P 



(a)(b)(c)(d)(e)(f)^a 



A (a) A W A ( C ) A W A (e) A (« 



(100) 
(101) 
(102) 



A a a ^ is constructed from an orthonormal set of vectors A^ (a = 0,1,2,3), which are parallel-propagated along a 
timclike geodesic Q as 



(103) 



Since the tangent vector of the timelike geodesic, u M = (d/dr)^, is parallel-propagated along Q, Aq should be equal 
to u M . Here, the geodesic equations are integrated to give (e.g., [15]) 



it* = ^[(J 2 + a 2 )Y,E - 2MarB], 
■ r " 2 - A 2 - A(r 2 +K), 

B 2 



{Hu e f =K-a 2 cos 2 - 



sin 2 9 



2 /)' 



(104) 
(105) 
(106) 
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1 



AE sin 2 



-2MrB + LE], 



where 



A = E{r 



2 ' a 2 ) 



aL, B = L — aE sin 



(107) 



(108) 



E = —Ut and L — u v are the specific energy and angular momentum for a particle of mass [x moving around a Kerr 
black hole. K = (L — aE) 2 + Ck is the so-called Carter constant with Ck a constant [25]. The first integral of the 
geodesic equations leads to 



,(°) 



A< 2) = VeV, 
A (3) = B 
V£sin0 



A 



VT,u r 
~7z' 



For other components of Ao°\ we follow Marck [15], and thus, we choose 



A< 0) 



a\/Eru r 
arA 



cos * 



aA 
v/AE 



sin*, 



cos * sin *, 



A ( 2) = _ PaBcosB ^ ^ _ ^ 

l(3) pVZau cos (3B 

A) ' = = cos * r= sin *, 

Vif VSsin6» 



(0) _ 

2 — 

(1) _ 



(3) 



A 



A (o) 



cos 6u r 

y/KA ' 

aA cos 
VXAE ' 

VXSsin^' 



'AT 
av / EY'u r 



A« = 



VXA 

orA 



Af) 



A 



(3) 



VifAE 
(3aB cos 
VAT] sin 

(3\/T,au e cos 



sin * H : cos *, 



/AS 

sm * H cos * 



VA 

sin* + /3V£u e cos*, 
0B 



where a and (3 are normalization constants defined by 



sin* 



y/Esm6 



cos*, 



K - a 2 cos 2 



/? = -■ 



(109) 

(110) 

(111) 
(112) 



(113) 
(114) 
(115) 
(116) 
(117) 
(118) 
(119) 

(120) 
(121) 
(122) 
(123) 
(124) 

(125) 



We note that the direction of the components 1, 2, and 3 [not (1), (2), and (3)] adopted by Marck are approximately 
equal to x, —z, and y of the Cartesian coordinates in the comoving frame. The time evolution of the rotation angle 
* is computed by 
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IV. FORMULATION FOR EQUILIBRIUM NEWTONIAN STARS IN EQUATORIAL CIRCULAR 

ORBITS IN THE BLACK HOLE TIDAL FIELD 



A. Basic equations 



In this section, we give a formulation for computing equilibrium states of a fluid star in circular orbits around 
equatorial plane of a Kerr spacetime. Here, we assume that the self-gravity of the star is described by the Newtonian 
gravity. In this case, the gravitational potential associated with the tidal potential can be linearly superposed [13]. 
Using this property, we write the hydrodynamic equation for the fluid body as 



dvj | v jdvj = dP d(cj) + (fttidai) 
dr dx J dx l dx l 



( dAj _ dAA _ dAj 
dx % dxi J Q T 



(127) 



where p is the mass density, v l the three-velocity {dx l /dr), P the pressure, and </> the Newtonian potential produced 
by a star which obeys the Poisson equation 



A<f> = 4irp. 



(128) 



Ai is a vector potential defined in Eq.(134), which is associated with the so-called gravitomagnetic force [26]. </>tidai 
denotes the tidal potential associated with the background Kerr spacetime, which is related to the metric computed 
in previous sections as 

^tidal = -^(.900 + 1) 

= ~goo,ijX l x j - -^gm,tjkx % x 3 x k - -^goo,ijkix l x : >x k x l + 0(x 5 ) 

= ijCyarV + \c ijk x i x :i x k + [C m + ±C {ij C kl) - 4B {klw B ij)n ]x i x j x k x l + 0(x 5 ), (129) 



where 



o?:oj . 



ijk 



C, 

Cijkl = Ro(i\0\j;kl) , 
Bijk = Rk(ij)0, 

2 

—j 

3 



0(*|0|j;fe)) 



o k X X 3 . 



(130) 
(131) 
(132) 
(133) 

(134) 



In the equations of motion (127), we include the lowest-order gravitomagnetic term associated with Ak, although 
it is a first post-Newtonian term and does not appear in Newtonian order from the point of view of post-Newtonian 
approximations [27]. The reason we add it is that the order of magnitude of this term is as large as that of the 
fourth-order terms in </> t id a i if the spin angular velocity of the fluid star is of order f2 as in the corotational velocity 
field (see a discussion in the final paragraph of this section). On the other hand, for the irrotational velocity field in 
which the spin of the fluid star is negligible in the frame of the Fermi normal coordinates, the magnitude of this term 
will be much smaller than the fourth-order term in </> t idai- In the following calculation, we neglect the gravitomagnetic 
terms for most of calculations, but to clarify the quantitative effect of this term, we perform a few computations 
including it. 

In this paper, we restrict out attention to circular orbits in the equatorial plane, i.e., 9 = ir/2, u r = 0, u e = 0, and 
C K = 0. Then [20] 



2 - 2Mr + aVMr 

Wr{r 2 - 2a\/Mr + a 2 ) 
rD ; 

D = \Jr 2 - 3Mr + 2aVMr, 




(135) 
(136) 
(137) 
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where r denotes the orbital radius. In this case, the evolution equation for W becomes 

[m , , T [M ,,„„, 

d7 = V r^' ' V ^ T - (138) 

For the equatorial circular orbits, the transformation matrix reduces to a simple form as 



A o a) = (\/i + ^r,o,o, 7 ), (is-)) 



B 2 B 



r 



A (o) = / _ ^ sin cos 0,-Jl + ^sintf ], ( 140) 



A^ o) = (0,0,1,0), (141) 



A 3 o) = ( — cos sin 0, \/ 1 + cos * ) , (142) 



r 

where B = L — aE = r{\[Mr — a)/D. Note that 1 + B 2 /r 2 may be written as A/D 2 . A known interesting property 
is that independent of the value of a, B/r = l/\/3 and A/D — 4/3 at ISCOs [13] at which r satisfies 

r 2 - 6Mr + 8M 1/2 ar 1/2 - 3a 2 = 0. (143) 

This implies that at the ISCO, r 3 dj and r 3 Bijk are independent of a (see below). 

To derive the tidal tensors dj, djk, Cijki, and Bijk, as a first step, it is better to calculate the components in a 
tetrad frame defined by 

A«=A<'\ (144) 

A^ Q) = A^ a) cos + A 3 a) sin tf , (145) 

A^=A 2 a) , (146) 

A 3 a) = -A[ a) sin * + A 3 a) cos * . (147) 

We refer to this frame as the tilde frame in the following. In the tilde frame, A^ is independent of 'J, but the 
coordinate basis of this frame is not parallel-transported along the timelike geodesic. In the second step, we should 
perform the coordinate transformation to the parallel-transported frame. 

The nonvanishing components of dj, djk, djkl, an d Bijk, which denote the tidal tensor in the tilde frame, are 

M ( r 2 + B 2 \ 

^1 = ^(1-3^^], (148) 

- M ( B 2 \ , s 

^2 = ^(^1 + 3^, (149) 

M 

dz = -r, (150) 



#131 = B311 = —B232 = -B322 = -^B 113 — -B223 — 2r^~ \ 1 + "r 2 "' (151) 

3MA 1 / 2 

dn = 7 [2Dr 2 -2aBr + W 2 D], (152) 
Dr' 

MA 1 / 2 

C122 - di2 = C221 = -^[3i> 2 - 8aBr + 7B 2 D], (153) 

Dr 1 

MA 1 / 2 

C133 = C 3 i 3 = C331 = ^^[3i> 2 + 2aBr + 2B 2 D], (154) 

Gnu = ^ [-8L>r 4 + \8DMr 3 + l&aBrA - \2{a 2 + B 2 )Dr 2 + 27 B 2 DMr - l9a 2 B 2 D] , (155) 
C1122 = C1212 = C1221 = C2112 = C2121 = C2211 

= ^9 [24Dr 4 - hlDMr 3 - l08aBrA + 51(a 2 + B 2 )Dr 2 - l09B 2 DMr + 102a 2 B 2 D] , (156) 
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Cll33 = C1313 — C1331 — C3113 — C3131 — C3311 

- M ■ [24Dr 6 - HlDMr 5 + 20aBr 3 A + 25(a 2 + B 2 )Dr 4 - 56B 2 DMr 3 



2L>r n 

+10aB 3 rA + 5B 4 Dr 2 + 2ha 2 B 2 Dr 2 - \hB 4 DMr + 10a 2 B 4 D] , (157) 
C2222 = [3i> 4 - 6DMr 3 - WaBrA + 7(a 2 + B 2 )Dr 2 - UB 2 DMr + 19a 2 B 2 D] , (158) 

C2233 = ^2323 = C2332 = C3223 = C3232 = C3322 

= -^-^ [-6Dr 6 + l2DMr 5 + WaBr 3 A - 10(a 2 + B 2 )Dr A + 2lB 2 DMr 3 

-10aB 3 rA - 5B 4 Dr 2 + 5a 2 B 2 Dr 2 + 15B 4 DMr - lOa 2 B 4 D] , (159) 

O 71 f 

C3333 = \-3Dr 4 + 6DMr 3 - 6aBrA - 3(a 2 + B 2 )Dr 2 + 7B 2 DMr - 6a 2 B 2 D] . (160) 
Dr J 

The expression for Cij agrees with that derived by Marck [15]. The nonvanishing components of Cy, Cijk, Cijki, and 
Bijk are derived by the coordinate transformation from x l to x % . 

In the tilde frame, the tidal potential up to the fourth order is written as 

^tidai = ^C ij x i x j + ^C ljk x t x 3 x k + ^ [C^ki + 4C {ij C k i) - AB {U \ n \B ij)n ] x i x j x k x l , (161) 



where 



r — x 1 cos* + x 3 sin*, (162) 



x 1 

Z2 — rJ- 



i z =x 2 , (163) 
x 3 = -x 1 sin* + x 3 cos*. (164) 

In the Newtonian limit r » M(> a), it is written as 

M 



A, - M 
Ptidal - ^3 



2{x l f + (x 2 ) 2 + {x 3 f -2(x') 2 + 3{(x 2 r + (x 3 ) 2 } 



M 
"8^ 



.r; 



2r 

) 4 + 3(x 2 ) 4 + 3(Z 3 ) 4 - 24{(i 1 ) 2 (i 2 ) 2 + (x 1 ) 2 ^ 3 ) 2 } + 6(x 2 ) 2 (x 3 ) 2 l . (165) 



This agrees with the expansion form of the Newtonian tidal potential from a point source of mass M at a distance r 
as 

= (166) 
v /(i 1 +r) 2 + (i 2 ) 2 + (i 3 )2 

The orders of the magnitude of the second-, third-, and fourth-order tidal forces in Eq. (165) are 0(MR/r 3 ), 
0(MR 2 /r 4 ), and 0(MR 3 /r 5 ), respectively. On the other hand, the order of magnitude of the gravitomagnetic force 
in the Fermi normal coordinates is 0(M 3 / 2 Rv/r 7 / 2 ) where v denotes the characteristic magnitude of v l . For the 
corotational velocity field, v — 0(M 1 ^ 2 R/r 3 ^ 2 ), and hence, the gravitomagnetic tidal force is of 0(M 2 R 2 /r 5 ) which 
is the same as the order of the fourth-order tidal potential for stars with R = 0(M). For close corotating orbits 
with r — O(M), it is as larger as the third-order term. For rapidly rotating stars with v < c, the gravitomagnetic 
tidal force is always larger than the third-order term. On the other hand, for the irrotational velocity field, it will be 
negligible. 



B. Hydrostatic equations for corotational and irrotational equilibria 

Now, we turn our attention to the hydrostatic equations. First, we consider the case in which the velocity field is 
corotational, and assume 

v l = [-{x 3 - x c sin(fh-)}, 0, {x 1 - x c cos(rJr)}] (167) 

where f2 = d* /dr =const. x c is a correction constant which is much smaller than the stellar radius and nonzero only 
when we take into account the third-order terms in i^tidai or the gravitomagnetic terms. For x c 7^ 0, the rotational 
axis deviates from the x 2 axis. Also, the center of mass of a fluid star is different from the origin slightly (see Sec. V). 
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In this velocity field, Eq. (127) is integrated to give 

n 2 



[{x 1 - XgY + {^Y] =h + <j) + </> tidal + </> mag + C, (168) 
where x g = 2x c , C is an integration constant, and 

h = J —, (169) 



MB I B 2 



taag = 2— 1 + — 



-(^ ) 3 +i l { (^) 2 _( i 3 ) 2 }+ - {( -l )2 _ ( -2 )2} 



(170) 



The first integral of the Euler equation is independent of r in the tilde frame. Equations (128) and (168) constitute 
the basic equations for the corotational binary. 

For the irrotational velocity field, we should set the three-components of the four-velocity as 

dip 

Ui = Vi+Ai = -^r, (171) 

where ip denotes the velocity potential which is time-independent in the tilde frame. Then, Eq. (127) is integrated to 
give 

dip 1 dip dip d ^Aj,n M7o\ 

~W 2^ WW = h + + feal _ W A + a (172) 



In the tilde frame, the first term is written as 



Also, we have a relation 



(hP_(hp__ (hp_(hP_ 



Thus, the first integral of the Euler equation is also independent of r in the tilde frame. 
ip is determined by solving the continuity equation rewritten as 

where A is the Laplacian in the tilde coordinates. Equations (128), (172), and (175) constitute the basic equations 
for irrotational binaries. 



V. ROCHE LIMIT IN EQUATORIAL CIRCULAR ORBITS 



To quantitatively illustrate the importance of the higher-order terms in the tidal potential </> t idai as well as to clarify 
the dependence of the tidal disruption limit of a star on the equations of state and on general relativistic effects of the 
black hole, we numerically compute corotating equilibria and determine the tidal disruption limit (Roche limit) as an 
extension of a previous work by Fishbone [13]. For some case such as binaries of a black hole and a neutron star, the 
irrotational velocity field is more realistic. However, we know that in the incompressible case, the tidal disruption 
limit depends weakly on the velocity profile as far as the spin angular velocity of the star is of order M 1 / 2 /r 3/ ' 2 [17]. 
We expect that this will be also the case for compressible stars, and hence, even in the assumption of the corotational 
velocity field, we can obtain an approximate result of the tidal disruption limit for the irrotational velocity field. 
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A. Setting and numerical method 



We adopt polytropic equations of state for the star as 

P = Kp 1+ ~ , and thus, h = n[n + l)p~ , (176) 

where k is the polytropic constant and n the polytropic index. In this paper, we choose n — 0.5, 1, and 1.5 to 
approximately model neutron stars or white dwarfs. 

The basic equations in this problem are Eqs. (128) and (168). Numerical solutions are obtained by iteratively 
solving these coupled equations. To achieve a convergence in the iteration, we rescale the coordinates as x % = pq l 
where p is a constant and q l dimensionless coordinates. Then, the basic equations are written in the form 

A g = 47rp, (177) 

\p\<t - q g f + (f) 2 } = h + P 2 4> + P 2 (0 tidal + <^ mag ) + c, (178) 

where A q is the Laplacian in the coordinates of q' 1 , 4> — p~ 2 <f>, </>tidai = P ^tidai, 0mag = P ~ 2 ^mag, and % = q x g . 
Thus, six free constants M, a, k, q g , p, and C are contained in the equations. In the following, M is fixed adopting 
the units c = G = M = 1. An equilibrium configuration is computed for fixed values of k, p c (the central density), a, 
and r(> R). Sequences of the equilibria are computed varying these parameters, k and p c determine the mass m and 
the radius R of a star, and so do the ratios Q — m/M and R/r. Note that for a given value of r, fl is determined to 
be (M/r 3 ) 1 / 2 in this problem. 

Three remained constants q g , p, and C are parameters determined at each step of the iteration from the following 
conditions: During the iteration, we require that p = p c and dp/dq 1 = at the origin (q 1 , q 2 , q 3 ) = (0,0,0). In 
addition, we fix the coordinates of the stellar surface along the i 1 axis (the axis connecting the origin and the center 
of a black hole) on the black hole side as (q s , 0, 0) where q s < 0. From these three conditions, the three free parameters 
are determined. 

If we include the third-order terms in the tidal potential, the center of mass of a star may be deviated from the 
origin, although for consistency, it should be located approximately there. To check that the deviation is much smaller 
than the stellar radius, we calculate the value of q 1 coordinate for a center of mass defined by 

(q 1 )^^ f d 3 qpq\ (179) 
m J 

where 

m = p 3 J d 3 qp. (180) 

We found that | (q 1 ) | is indeed much smaller than the stellar radius (e.g., Kg 1 )! ~ 0.005|q s | at ISCOs for a = 0-0. 9M 
and the value is smaller for smaller orbital radii). 

The Poisson equation (177) is solved in the Cartesian coordinates with the uniform grid of size (2A+1, N+l, 2N+1) 
for (q 1 , q 2 , q 3 ) which covers the region with — L < q 1 < L, < q 2 < L, and — L < q 3 < L (the reflection symmetry 
with respect to the q 2 = plane is assumed) . The numerical method is essentially the same as that in [28] : We make 
the second-order finite-differencing equation and solve it in a preconditioned conjugate gradient method. Typically, 
N and the grid spacing A are set to be 50 and |(7 s |/40. To check the convergence, we varied the values of (N, \q s \/A) 
as (60,48), (50,30), and (40,32). It is found that the numerical results converge at the second order and the error is 
within 0.1% with (50,40). 

Following Fishbone [13], we often refer to a nondimensional parameter defined by 



irp c np c r A 



The order of magnitude of this parameter is 



(182) 

mr m/R z 

Thus, it denotes the ratio of the tidal force by a black hole to the self-gravity of a star. In particular, we focus on the 
value of C at the Roche limit, which is denoted by £crit in the following. £ cr j t is a function of r/M and depends on 
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the equations of state and the black hole spin. We can also compute the minimum allowed mass of a star that can 
escape tidal disruption for a given radius. The critical mass ratio associated with such minimum mass is defined by 
Qcrit = (m/M) minimum- Qcrit is also a function of v and depends on the equations of state and the black hole spin. A 
star will be tidally disrupted for Q < Q cr i t . 

Numerical computations for determining ( cr - lt and Q C rit are performed approximately fixing the value of an averaged 
stellar radius for a given set of r, a, and n. The stellar radius of a spherical polytrope Rq is written by [29] 



Ro 



(n + l)np, 



(l-n)/n 



4ir 



1/2 / , X 1/2 



a = (^) (183) 



where h c = n{n + l)p c and £x denotes the Lane-Emden coordinate at the stellar surface, which is 2.75270, 7r, and 
3.65375 for n = 0.5, 1, and 1.5, respectively. Note that for n = 1, fixing the value of n is equivalent to fixing the value 
of Ro. For other values of n, k varies along a sequence of a fixed value of Ro. 

To determine the Roche limit for a given set of r, a, and n, we compute a sequence of solutions by varying p c from 
a large value to a small value until an inner edge of the star forms a cusp on the black hole side. A configuration 
with such a cusp can be identified as the Roche limit, i.e., the self-gravity of the star is small enough to form a saddle 
point of the total gravitational potential at a stellar surface. Specifically, the Roche limit is determined monitoring 
the following quantity at (q 1 , q 2 7 q 3 ) — (q s ,0, 0): 

H = n (q s - q g ) — . (184) 

Here, H is proportional to dh/dq 1 . Thus, the value of H for q s < is positive for stable stars, and becomes zero for 
the marginally stable configuration against tidal disruption. 

In this paper, we are in particular interested in making approximate models of binaries composed of a stellar-mass 
black hole of mass M <~ 3-30M Q and a neutron star of mass m ~ 1.4M Q . It is appropriate to assume that the radius 
of neutron stars is between 10 and 15 km. Then, Ro in the present units should be between ~ 3M and ~ 0.2M. In 
the following calculation, we adopt the values of Ro/M as 0.5, 1, and 2. For consistency in the framework of this 
paper, (i) Rq should be much smaller than r, (ii) m should be much smaller than M, and (iii) r should be larger than 
the orbital radius of an ISCO (hereafter nsco) around the black hole. Thus, the computations are performed only for 
m < M and for r > nsco ■ ^ would not be quantitatively appropriate to model a neutron star by Newtonian gravity 
since neutron stars are compact with m/Ro ~ 0.2 and thus general relativistic objects. The following numerical 
results could contain a systematic error of magnitude <~ m/i?o- However, the quantitative importance of the higher- 
order terms as well as the general relativistic effects in the tidal potential can be clarified even if we neglect the 
general relativistic corrections of neutron stars. Also, the present study will be useful for qualitatively clarifying the 
dependence of the Roche limit on the equations of state. 



B. Tidal potential 

In Fig. 1, we display the tidal potential ^tidai along x , x 2 , and x 3 axes. Figures l(a)-(c) show ^tidai for 
(r/M,a/M) = (6,0), (2,1), and (6,1), respectively. Figure 1(d) is 0tidai m the fourth-order approximation for 
a = M and r/M — 1.2, 1.5, 2, and 3. To clarify the convergence with increasing the order, ^tidai in the second- 
(dotted-dashed curves), third- (dashed curves), and fourth-order (solid curves) approximations are shown together in 
Figs. l(a)-(c). </>tidai in the second- and third-order approximations are identical along x 2 and x 3 axes. Thus, only 
the second-order results are presented. At the second order, 0tidai is symmetric with respect to the origin along all 
the axial directions. The third-order terms induce an asymmetry in the x 1 direction. 

The difference in the magnitude of </>tidai is about 20-30% between the second- and fourth-order approximations for 
\x l \ ~ Ro ~ M near the ISCOs. For \x l \ < R - M, the difference in the magnitude of the third- and fourth-order 
tidal potentials is < 1%, and hence, the convergence is approximately achieved. </> t idai in the third- and fourth-order 
approximations are larger than that in the second-order approximation in the black hole side (x 1 < 0). This indicates 
that a star becomes prone to tidal disruption in the higher-order approximations. 

Comparing Figs. 1(a) and (c) for which the results with the same value of r are shown, it is found that the spin 
effect reduces the magnitude of </>tidai along x 1 and x 2 axes. This illustrates the property that for the larger value 
of a/M, the tidal effect is weaker for a given value of r/M. Figure 1(d) shows that the magnitude of r 3 t idai in the 
black hole side is smaller for the smaller values of r. This suggests that a star is less prone to the tidal disruption for 
an extremely high value of a w M near the ISCOs. Such a special feature is not outstanding for other values of a. 



18 



C. Roche limits 



In Fig. 2, we show Ccrit and /i cr ; t = Q C rit{r/Ro) 3 as functions of r/M (a) for n = 1, a = 0, and Ro/M = 0, 0.5, 1, 
and 2, and (b) for n = 1, a = 0.9M, and Ro/M = 0, 0.5 and 1. Stars with C < Ccrit and fj, > /i cr it for a given value of 
r/M are stable against tidal disruption. The computations for R ^ were performed taking into account the tidal 
potential up to fourth order. Note that the values of Ccrit and /x cr it m the third- and fourth-order approximations with 
i?o = are equal to those in the second-order approximation with an arbitrary value of Ro- Thus, 11 Ro = 0" implies 
that the computations were performed in the second-order tidal approximation. With increasing Ro/M, Ccrit (/i C rit) 
at a given orbital radius decreases (increases). This illustrates that the tidal force in the fourth-order approximation 
is stronger than in the second-order one. 

As in previous works [13,17] carried out by the second-order tidal approximation, Ccrit (Merit) decreases (increases) 
with decreasing r/M. This implies that with the decrease of r/M, the tidal force is enhanced. At r = nsco, Ccrit 
and /z cr ;t become minimum and maximum for given values of a and R . We denote them by Ccrit:min and ^ C rit:max- 
On the other hand, for r — > oo, Ccrit and ^ cr it become maximum and minimum for given values of a and Ro, and they 
agree with the values in the Newtonian limit [18]. For given values of a and Ro, a star at an ISCO has to be tidally 
disrupted whenever C > Ccrit > Ccrit:min or fj, < ^ crit < /i crit:mm . In other words, for C < Ccritimin and /x > Mcrit:max, the 
star in a circular orbit is not tidally disrupted outside the ISCOs. 

In Figs. 2(c) and (d), Ccrit /Ccrit :2nd - 1 and Merit /Mcrit:2nd — 1 as functions of M/r are shown for the same parameters 
as in Figs. 2(a) and (b), respectively. Here, Ccrit:2nd and Mcrit:2nd denote Ccrit and ^ crr t determined in the second-order 
tidal approximation. We note that in the second-order tidal approximation, Ccrit and ^ cr it are independent of Ro. 
It is found that for M/r(< 0.1), Ccrit /Ccrit :2nd — 1 and Mcrit/Mcrit:2nd — 1 are approximately linear functions of M/r. 
Thus, for the large values of r/M, the third-order tidal potential dominantly modifies the values of Ccrit and /z cr it- 
The numerical results show that the following relations approximately hold: 



Here, the coefficients of the correction factors C\ and C2 are w 0.95 and w 0.80 irrespective of the value of a for n = 1. 
The error is within ~ 0.05 for both coefficients. For M/r > 0.1, the fourth-order tidal potential becomes important. 
However, Eqs. (185) and (186) approximately hold even for r > 6M. For r < 6M, the correction factors in these 
fitting formulae give overestimated values. 

To clarify the importance of the higher-order terms of Ro/r in the tidal potential, in Fig. 3, we show Ccrit and ^ crr t 
as functions of r/M (a) for n = 1, a = 0, and Ro = M, and (b) for n = 1, a = 0.9M, and Ro = M with the second-, 
third-, and fourth-order tidal potentials (dotted, dashed, and solid curves, respectively). With the higher-order tidal 
potential, the value of Ccrit decreases. Namely, the minimum allowed value of p c increases. This also illustrates that 
in the second-order approximations, the tidal force is underestimated. For r/M{— r/i? ) = 6 (10), the values of Ccrit 
in the third- and fourth-order approximations are about 13 (8) and 15 (10)% smaller than that in the second-order 
one. This shows that the standard tidal approximation taken up to the second-order term provides a result with the 
error of ~ 1Q{Rq/M)% for close orbits with Ro/r > 0.1. However, the convergence appears to be very good if the 
third- and fourth-order terms are included for R < M. In other words, the fourth-order term plays a quantitatively 
minor role. All these results agree qualitative with a Newtonian analysis [18]. 



n 


Ccrit:2nd(ISCO) 


Mcrit:2nd(ISCO) 


Ccrit ( r -> 00) 


Merit (r — ► 00) 


M C (M Q ) 


Ci 


C 2 





0.0664 


20.1 


0.0901 


14.8 


4.54 






0.5 


0.0480 


19.4 


0.0651 


14.5 


4.45 


0.7 





1.0 


0.0303 


14.9 


0.0411 


11.1 


3.90 


0.95 


0.80 


1.5 


0.0171 


13.7 


0.0232 


10.1 


3.74 


0.97 


0.97 



TABLE II. Values of Ccrit:2nd, Mcrit:2nd, M c , C\ , and C2 for each value of n. The values for n — were determined by 
Fishbone [13]. Ccrit:2nd and /i cr it:2nd do not depend on the value of a, and Ci and C2 depend very weakly on it. The value of 
M c shown here is that for a = and r = nsco — 

6M, and for a ^ 0, M c (a) = M c (a = 0)(6M/nsco) 3/2 . For the values of Ci 
and Ci, the numerical error is within 5% for n = 1 and 1.5, and ~ 10% for n = 0.5. 




(185) 



(186) 
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In Fig. 4(a), we show £ cr i t and /i cr ;t as functions of r/M for n = 1, i?o = M, and a/M = 0.9, 0.8, 0.3, 0, —0.3, and 
— 1 in the fourth-order tidal approximation. For comparison, the results in the second-order approximation are shown 
in Fig. 4(b). These figures clarify effects of the black hole spin on the Roche limit. For each curve of a given value 
of a, the minimum value of r is equal to nsco- ^ n tne second-order tidal approximation, the values of ( a a and yu crit 
at the ISCO are independent of a [13] and are about 0.0303 and 14.9, respectively (cf. Table II). In the third- and 
fourth-order approximations, on the other hand, these values depend slightly on the value of a. As mentioned before, 
the value of £ C rit (/i C rit) in the fourth-order tidal approximation is smaller (larger) than that of Ccrit:2nd (Mcrit:2nd)- The 
magnitude of the difference between £ cr it and (crit:2nd is slightly smaller for a large value of a > 0.9M and a < 0. The 
effects of the third- and fourth-order terms in the tidal potential are strongest for 0.3 < a/M < 0.8. 

In Figs. 5(a)-(d), the density contour curves in the equatorial plane (x 1 -^ 3 plane) at the Roche limit are displayed. 
Figures 5(a)-(c) are those in the fourth-, third-, and second-order tidal approximations for r = 6M, i? = M, a = 0, 
and n = 1, and Fig. 5(d) is in the fourth-order approximation for r = 2AM, Rq — M, a = 0.9M, and n = 1. In the 
second-order approximation, the density is symmetric with respect to the x 1 = plane, and hence, the asymmetry 
is induced by the third-order term. At the tidal disruption, x s = pq s w 1.6, 1.55, and 1.5i?o in the second-, third-, 
and fourth-order approximations. On the other hand, the length of the minor axes is ~ Ro independent of the order 
of the approximation. Thus, in the higher-order approximations, the ellipticity at the Roche limit is slightly smaller. 
Comparing Figs. 5(a) and (d), it is found that the density configurations near ISCOs for different values of a are very 
similar irrespective of r. 

To clarify the dependence of the Roche limit on the equations of state, we show £ crit and ^ cr i t for n = 0.5, 1, and 1.5 
and (a) for a — and (b) for a — 0.9M in Fig. 6. The solid and dotted curves for each value of n denote the results in 
the fourth-order tidal approximation with Rq = M and in the second-order tidal approximation, respectively. Recall 
that for ( > ( cr [ t or \i < jU cr i t , a star is stable against tidal disruption. This implies that for given values of moss and 
radius (Rq < M), a star with softer equations of state is stronger against tidal disruption than that with the stiffer 
one. On the other hand, £ cr it is smaller for softer equations of state. This implies that for given values of central 
density and radius (Ro ;$ M), a star with stiffer equations of state is stronger against tidal disruption * 

From the analysis of Ccrit/Ccrit:2nd — 1 and ^ cr i t /^ C rit:2nd — 1 as functions of M/r, the coefficients C\ and C2 are 
computed for a = and 0.9M. We find that they depend very weakly on a and C\ <~ 0.7 and C2 ~ for n = 0.5, and 
C\ w 0.97 and C2 ~ 0.97 for n = 1.5. For n = 0.5, C\ and C2 are not determined very accurately, and hence, the error 
would be ~ ±0.1. This is because for very stiff equations of state, the density steeply decreases near the stellar surface 
and it is not easy to accurately apply the condition for determining the Roche limit (cf. Eq. (184)). Nevertheless, 
Ci and Ci for n = 0.5 are much smaller than those for n = 1 and 1.5, and therefore, we can conclude that C\ and 
C2 are larger for softer equations of state. This implies that the third- and fourth-order terms in the tidal potential 
play a more important role in softer equations of state. For n = 0.5, C2 is approximately zero within the numerical 
error. This can be explained as follows. The central density at the Roche limit in the fourth-order approximation is 
always larger than that in the second-order one. On the other hand, the volume of a star at the Roche limit in the 
fourth-order approximation is smaller than in the second-order one. These two effects seem to approximately cancel 
for n — 0.5. 

Computations were also performed including the gravitomagnetic terms associated with Ai for n = 1 , Rq /M = 0.5,1, 
and a/M = 0, 0.9. In Fig. 7, we display C cr it and /i cr it as functions of r/M. This figure shows that the gravitomagnetic 
terms decrease Ccrit and increase yu cr i t for a given value of r, respectively. This implies that the gravitomagnetic tidal 
force reduces the magnitude of the tidal force. This result is reasonable since the gravitomagnetic force in the Fermi 
normal coordinate is mainly associated with the coupling between the star's spin and orbital motion, and the spin- 
orbit coupling force has a repulsive nature in the case that their axes are parallel [30] . The order of magnitude of this 
term is smaller by the order of r~ 2 as described in the end of Sec. IV A. Thus, it plays an important role only for 
very small orbital radii r < 5M as in the fourth-order term of <^tidai- The numerical results show that at r = 6M, 
Ccrit and ^ cr it change only by ~ A(R a /M)% for a = and by ~ 3(Ro/M)% for a = 0.9M. However, for rapidly 
spinning black holes, this term can significantly modify the Roche limit near the ISCOs. Indeed, for a — 0.9M with 
t ~ nsco ~ 2.32M, their values are changed by 20-30% for Rq w M. 



^Figure 6 suggests that these statements may not be correct for a very large value of Ro S> M. However, for such a large 
value of _Ro, the tidal approximation is not applied, and thus, this point is not clear. 
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D. Application to neutron star-black hole binaries 



Now, we apply the results obtained in the previous subsections to binary systems of a neutron star and a black 
hole. Using the result for /z cr it = /U C rit:2nd(l + C2-R0A) = (rn/M) CI it{r/Ro) 3 , the condition for the mass of a black 
hole that can tidally disrupt a neutron star outside the ISCO is derived as 



where M c {a) is a function of a. M c is estimated for r = nsco, and thus, it denotes the maximum mass of a black hole 
that can tidally disrupt neutron stars of a given set of m = 1.4M©, R = 10 km, and nsco/M. The dependence of M c 
on a only results from the dependence of nsco /M on a since /x cr it:2nd at the ISCO is independent of a (/i C rit:2nd = 19.4, 
14.9, and 13.7 at the ISCOs for n = 0.5, 1, and 1.5; cf. Table II). 

For a = 0, M c (a = 0) w 4.45M , 3.90M©, and 3.74M© for n = 0.5, 1, and 1.5, respectively. Note that M c (a = 
0) w 4.54M Q and /i cr it:2nd = 20.1 at ISCOs for n = [13,17], and thus, they are close to the values for n = 0.5. M c (a) 
for a 7^ can be computed by M c (a = 0)(6M/nsco) 3 ^ 2 - For a rapidly rotating black hole with a — 0.99M (0.9M), 
nsco/M w 1.4545 (2.3209), and hence, M c w 37.3M Q , 32.7M , and 31. 3M© (18.5M Q , 16.2M , and 15.5M ) for 
n = 0.5, 1, and 1.5, respectively. Thus, for the large value of a/M = 0.9-1, neutron stars of i?o ~ 10 km are tidally 
disrupted for a wide mass range of the black hole with M < 15-40M Q . If the typical radius of neutron stars is i?o = 15 
km, the value of M c increases by a factor of w 1.84, and M c becomes s=s 30-75Mq for a/M = 0.9-1. 

In the second-order tidal approximation, M c for given values of Rq and m is smaller for the larger value of n. It 
is interesting to point out that difference between M c for n = 0.5 and for n = 1 is fairly large ~ 14% although the 
differences for n — 1 and n = 1.5 and for n = and n = 0.5 are only ~ 4% and ~ 2%, respectively. This suggests 
that the critical mass ratio <2 cr it for orbits close to the ISCO depends sensitively on the equations of state in a stiff 
region of n = 0.5-1. 

For n = 1 and 1.5, the higher-order terms in the tidal potential increases the critical mass for the tidal disruption. 
For parameters M = 4.5M Q and i?o = 10 km, we obtain Rq rj 1.5M. In this case, the critical mass with a = and 
nsco = 6M increases by a factor of « 1 + C2/8. For M = 15Mq and Rq = 10 km, Rq w 4M/9. In this case, the 
critical mass with a — 0.9M and nsco ~ 2.32M increases by a factor of rj I + C2/IO. Thus, the critical mass for orbits 
close to the ISCO increases by ~ 5% beyond the value of M c for n = 1 and 1.5 due to the effect of the higher-order 
tidal potential, and hence, the dependence of the critical mass on n would be weaker in reality. Nevertheless, the 
critical mass near the ISCO still depends fairly strongly on n for n = 0.5-1 which are plausible polytropic indices for 
neutron stars [29]. In contrast, £ cr i t is smaller for softer equations of state, and the higher-order terms in the tidal 
potential make this feature stronger. 

To clarify its dependence on the orbital radius, in Fig. 8, we display the critical mass of a black hole as a function 
of r/M. Here, the critical mass in the fourth-order approximation is defined by 



For M < M cr n at a given orbital radius r/M, the neutron star is tidally disrupted. In Fig. 8, we show M cr ; t for a = 
and 0.9M, for R = 10 and 15 km, and for n — 0.5, 1, and 1.5. It is found that if we assume that the mass of black 
holes is larger than 3M Q , the tidal disruption can happen only for close orbits with r/M < 7 for R = 10 km and for 
r/M < 10 for Rq — 15 km. Based on an observational results for black hole binaries in our galaxy and in the LMC, 
typical mass of black holes is in the range between 6 and 8M Q [33] . This suggests that the tidal disruption of neutron 
stars may happen frequently only if (i) the radius of neutron stars is fairly large as ~ 15 km or (ii) the typical spin 
parameter of black holes is fairly large as a/M > 0.6. 

As mentioned above, the values of M cr ; t are very close each other for n = 1 and 1.5 irrespective of the orbital radius 
r/M, and difference in the value of M cr j t for n = 0.5 and 1 is remarkable for r ~ nsco- It is interesting to note that the 
value of Merit depends very weakly on the equations of state for orbits not very close to the ISCO (or in other words for 
M cr ;t ^ 4M Q ). The reason is that with increasing r/M, the value of M cr it steeply decreases as M„i t /R oc (M/r) 3 / 2 , 
and hence, for a given value of Rq, the importance of the correction due to the higher-order terms associated with the 
term Ro/r — (Ro/M)(M/r) in determining M cr ;t is enhanced for n = 1 and 1.5 § . Nevertheless, M cr ;t for small values 



§ The tidal approximation cannot be used for M cr i t < 2M because M cr i t and r should be much larger than m and 7?o, 
respectively. Thus, one should focus only on the results for M > M cr i t in Fig. 8. 




(187) 




(188) 
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of r/M < 5 (or for large black hole mass with M > 4M ) still depends on the equations of state. This suggests that 
by determining the tidal disruption limit of binaries of a neutron star and a massive and rapidly spinning black hole 
in an observation, the equations of state of neutron stars may be constrained. 

An observation of binaries composed of a stellar-mass black hole and a neutron star will be possible in the near 
future using laser interferometric gravitational wave detectors such as LIGO [31]. For the orbital separation r > 10M, 
the binary adiabatically evolves due to emission of gravitational waves. In such inspiral phase, the masses of the black 
hole and the neutron star as well as the black hole spin will be determined from the chirp signal of gravitational waves 
by using matched filtering techniques [32]. If the mass of the black hole is smaller than the maximum value of M crrt , 
the tidal disruption of the neutron star will happen near an ISCO. It is reasonable to expect that the chirp signal 
of gravitational waves quickly shutdowns at the tidal disruption. This suggests that from the signal of gravitational 
waves, the radius and the orbital frequency at the tidal disruption may be identified. With increasing observational 
samples, M cr ;t (or Q C rit) may be determined. If that becomes possible, the equations of state of neutron stars will be 
constrained. 

Unfortunately, the present formulation is not yet appropriate for an accurate determination of the values of M cr it 
and Qcrit for a binary of a neutron star and a black hole, since m and Ro are not much smaller than M and r. Although 
it is reasonable to expect that dependence of M cr ;t and Q C rit on the equations of state is qualitatively unchanged even 
in a more accurate computation, the values would be systematically modified by 10-20% (see Appendix A). Thus, 
M cr it and Qcrit should be determined more accurately in terms of fully general relativistic computations in the future. 

VI. SUMMARY AND DISCUSSION 

As an extension of previous works by Fishbone [13] and Marck [15], we have derived the tidal potential induced 
by a black hole up to the fourth order in R/r in the Fermi normal coordinate system using the method developed 
by Manasse and Misner [19]. The new tidal potential is incorporated into the Newtonian equations of motion for 
a star orbiting the black hole. Using the new formulation, we determined the tidal disruption limit (Roche limit) 
for corotating Newtonian stars in equatorial circular orbits around a black hole. It is found that the third- and 
fourth-order terms in the tidal potential always amplify the tidal force and modify the Roche limit for close orbits. In 
particular, the third-order term plays a quantitatively important role for orbits with R/r > 0.1. To this time, tidal 
problems for a star orbiting a black hole have been widely studied taking into account only the second-order terms 
of the black hole tidal field (e.g., [13,9,16,17]). The present results illustrate that for close orbits with R/r > 0.1, the 
second-order approximation might not provide quantitatively accurate results. 

For a specific illustration of the importance of the higher-order terms in the tidal potential as well as the dependence 
of the Roche limit on the equations of state of the star, numerical computations are performed for plausible parameters 
of binaries of a black hole and a neutron star. Since neutron stars are general relativistic objects and the mass ratio 
Q between two stars arc not very small, the present framework might not be yet appropriate for such studies. 
However, it is still possible to extract many qualitative features on the tidal disruption limit. The following is the 
summary of the results: (i) general relativistic corrections amplify the tidal potential with the decrease of the orbital 
separation irrespective of the order of the tidal approximation; (ii) as found by Fishbone [13], in the second-order 
tidal approximation, the minimum value of £ cr it and the maximum value of /x cr ;t are independent of a. In the third- 
and fourth-order approximations, they depend on a; (iii) because of the third- and fourth-order terms in the tidal 
potential, the critical mass ratio Q cr it for the tidal disruption is changed by <~ 10-15% for close orbits near ISCOs 
with Ro ~ M; (iv) with the increase of the value of spin parameter a, the magnitude of the tidal potential decreases 
for a given value of r/M, and as a result, £crit increases and /z cr it decreases. This property is independent of the order 
of the tidal approximation; (v) for given values of central density and stellar radius (Rq < M) , a neutron star with 
stiffer equations of state is stronger against tidal disruption. On the other hand, for given values of stellar mass and 
radius {Rq < M), a neutron star with softer equations of state is stronger against tidal disruption. The maximum 
value of /i. cr it depends sensitively on stiffness of the equations of state for n — 0.5-1. 

As mentioned above, the present formulation is not yet appropriate for an accurate determination of the tidal 
disruption limit of a neutron star by stellar-mass black holes, since the effects of mass, spin, and multipole moments 
of the neutron star to the orbital motion [34,35] as well as its general relativistic self-gravity are neglected in the 
analysis. In Appendix A, we estimate the order of magnitude of the error associated with such neglected effects. 
To more accurately determine the tidal disruption limit, fully general relativistic computation is obviously necessary. 
We expect that our result presented in this paper could be a guideline for the future computation. In particular, we 
note that for the case that the mass and radius of a neutron star are much smaller than the black hole mass and 
orbital radius, respectively, the present formulation provides an accurate result for the Roche limit. The fully general 
relativistic results computed in the future should be compared with our numerical results to check the accuracy. 
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So far, we have focused on binaries of a black hole and a neutron star. The results with n = 1.5 may be used 
for determining the Roche limit of a white dwarf near the ISCO of a massive black hole. For typical values of mass 
<~ 0.7 'Mq and radius <~ 10 4 km of a white dwarf with n ~ 1.5 [29,36], Eq. (187) is written as 

/ n \ 3/2 / v. -1/2 / \ -3/2 / s 1/2 

M<l.axl& M J*) -£-) J-) £*) . (189) 



10 4 k m y \0.7Mq) \6M J \13.7 

In this case, Ro/r <C 1, and therefore, the third- and fourth-order terms in the tidal potential are not important and 
can be neglected. Equation (189) implies that an intermediate-mass black hole of M ~ 1O 3 M0 will tidally disrupt a 
typical white dwarf in a circular orbit at r ~ 160M. To tidally disrupt a typical white dwarf in a circular orbit, a 
supermassive black hole of M ~ 1O 6 M has to be rapidly rotating with a > 0.95M for which nsco ^5 2M. 

White dwarfs orbiting a supermassive black hole in galactic centers are likely to have highly elliptic orbits with 
E :=y 1 or parabolic orbits with E = 1 (e.g., [37,38]). Even in this case, the values of M cr ; t , which are determined 
for the circular orbits with r w nscOj may be used for determining the upper mass of a black hole, M max , for which 
a white dwarf with E w 1 is tidally disrupted. The reason is that the tidal disruption limit will be determined by 
the tidal force at the periastron radius and thus M max will be determined by the radius of the marginally bound 
orbits. For highly elliptic or parabolic orbits in general relativity, the marginally bound orbits become the so-called 
zoom- whirl orbits [39], which have a nearly circular trajectory near the periastron with the orbital radius [29] 



r = r mb = 2M - a + 2M yjl - a/M. (190) 

Therefore, the analysis in assumption of the circular orbits is approximately applicable. Here, the circular orbits with 
r = r m b is unstable. Nevertheless, we can estimate the Roche limit by the analysis presented in Sec V mathematically. 

The result is that ^ C rit:2nd ~ 20.7 for n = 1.5 irrespective of value of a. Here, we note that for such circular orbits 
with E = 1 and r = r m b, B/r = 1 independent of a, and so is ^ C rit:2nd- Replacing nsco to r m b and adopting the new 
value of /i C rit:2nd, we can approximately estimate M max of the black hole for the tidal disruption of a white dwarf in 
highly elliptic or parabolic orbits as 

/ D \ 3/2 / \ -1/2 / v. -3/2 

M maxRj 3.77xlO^M (*) r^") ( r J±) . (191) 



^ 10 4 km y \O.7M J V4M, 

Therefore, a black hole of M < 3.8 x 10 5 Af Q and a = (M < 3 x 10 6 M Q and a rts M) can tidally disrupt white dwarfs 
of m w 0.7M Q and R w 10 4 km. 

In this paper, we have only studied the tidal disruption limit for a star in equatorial circular orbits. The formulation 
derived in this paper can be used for the hydrodynamic tidal problem of an ordinary star or a white dwarf in parabolic 
orbits around a supermassive black hole [16]. In Appendix B, we write the tidal potential for a parabolic orbit in the 
equatorial plane, which may be used for an extension of the works presented in [16]. To confirm the prediction (191) 
for the tidal disruption of white dwarfs by a supermassive black hole, we plan to perform numerical simulations. 
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APPENDIX A: ORDER OF MAGNITUDE OF THE ERROR 



In the analysis for the tidal disruption limit of a star by black holes in terms of the present tidal approximation, we 
neglect the effects associated with the mass, spin, and multipole moments of the companion star to its orbital motion. 
Here, we estimate the order of magnitude of the error due to neglecting these effects. 

Neglecting the mass of the star results in the error of the orbital angular velocity by a factor of m/M. This error 
is included in the centrifugal force associated with the corotating velocity field (cf. Eq. (168)). Since the order of 
magnitude of the centrifugal force is nearly identical with the tidal force, the values of Ccrit and /i cr ;t which characterize 
the tidal disruption limit would be modified by a factor < m/2M — 1-20% for m = 1.4M Q and M = M crit ~ 4-50M Q . 
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An equation of motion for an extended body was considered in [35] ignoring the self-gravity of the body. If we take 
into account the spin and quadrupole moment of the body, we have the equations of motion 

£p a = \v p S^R aPlS + \j^V a R 0l5 z , (Al) 

where p a is the momentum vector, s an afhne parameter (see Eq.171 in [35]), S al} the spin tensor, and J a P~< & the 
quadrupole mass distribution of the extended body. In the post- Newtonian approximation, two terms in the right-hand 
side of Eq. (Al) are estimated as 

SFi s J«<V„fl w , ~ . (A3) 

where we assume that the orbital angular velocity f2 is equal to the spin angular velocity observed by a comoving 
flame, and that the quadrupole moment is induced by the black hole tidal field. The ratio of these terms to the 
Newtonian force is 

M _ (A4) 



Mm/r 2 r 3 
S\F%\ _ MR 
Mm/r 2 r 5 



(A5) 



Thus, these corrections modify the orbital angular velocity by a factor of <~ MR 2 /r 3 and ~ MR 4 /r 5 , respectively. 
They are likely to be much smaller than m/M for neutron star binaries. We note that for the irrotational velocity 
field, the correction due to the spin is approximately zero. 

With the modification of the orbital motion, the location of the ISCO will be modified. By this effect, £ C rit:min and 
/^critimax will be also modified. 



APPENDIX B: TIDAL TENSORS FOR EQUATORIAL PARABOLIC ORBITS 

In Sec. IV, we derive the formulation of the tidal problem for a star in equatorial circular orbits. Another interesting 
case is a parabolic encounter of an ordinary star or a white dwarf with a supermassive black hole. Here, we write the 
tidal potential for equatorial parabolic orbits. 

For equatorial parabolic orbits with E = 1 and L > L cr i t where L crlt is a critical value which depends on a and 
L„it = 4M for a — 0, the first integrals of the geodesic equations are written as 



(r z + a z )r z - 2Mar£ , (Bl) 
2MI + Lr 



Ar 2 

u r = ±V 1 , (B2) 



Ar 

and u 8 = 0. Here, 



V,= 



(B3) 



and £ = L — a. Then, the nonzero components of the tidal tensor in the tilde frame are 
M ( r 2 +( 2 \ 

Cii = ^r U- 3 — - (B5) 



M ( l 2 \ 

^22 = ^(1 + 3^), (B6) 
C33 - ^, (B7) 
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Components in the Fermi normal coordinates are calculated by operating rotational matrices associated with W as in 
Sec. IV. Here, evolution equation ol the rotation angle is written as 



d* L 
dr ~ r 2 + £ 2 ' 



(B24) 



If computations are performed in the tilde frame, *S> is not necessary for computing the tidal potential. Instead, d^/dr 
and d 2 ^ /dr 2 appear in computing the inertial forces in the equations of motion. 
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(c) (d) 

FIG. 1. The profiles of the tidal potential 0tidai along x 1 , x 2 , and x 3 axes (a) for r = 6M and a = 0, (b) for r = 2M and 
a — M, (c) for r = 6M and a = M, and (d) for r/M — 1.2, 1.5, 2, and 3 and a = M. The solid, dashed, and dotted-dashed 
curves in panels (a)-(c) denote 0tidai in the fourth-, third-, and second-order approximations, respectively. In panel (d), the 
tidal potential in the fourth-order approximation along x 1 and x 3 axes are shown. The units of G = M = 1 are adopted in 
these figures. 
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(c) (d) 
FIG. 2. (a) £ cr it and /u cr it as functions of r for Ro/M — 0, 0.5, 1, and 2, and for a = 0. (b) the same as (a) but for Ro/M — 0, 
0.5, and 1 and for a = 0.9M. (c) Ccrit/Ccrit:2nd — 1 as a function of M/r for Ro/M =0.5 (dashed curve) and 1 (solid curve), 
and for a — 0. (d) the same as (c) but for a = 0.9M. Here, n = 1. Ccrit:2nd is Ccrit with Ro = and equal to that in the 
second-order tidal approximation. We note that a star with £ > (ci-it ( or A* < Merit) for a given value of r/M is unstable against 
tidal disruption. 
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(a) (b) 

FIG. 3. Ccrit and (j, CT n as functions of r (a) for a = and (b) for a — 0.9M. For both cases, n = 1 and 7?o = M. The solid, 
dashed, and dotted-dashed curves denote the results in the fourth-, third-, and second-order tidal approximations, respectively. 




(a) (b) 

FIG. 4. £ C rit and /i cr it as functions of r for n — 1 and various values of a/M. (a) The dotted, dotted-dashed, dashed, 
solid, dotted-long-dashed, and long-dashed curves denote the results in the fourth-order tidal approximation with Ro — M for 
a/M = 0.9, 0.8, 0.3, 0, —0.3, and — 1, respectively, (b) The dotted, dotted-dashed, solid, and long-dashed curves denote the 
results in the second-order tidal approximation for a/M — 0.9, 0.8, 0, and —1, respectively. Note that in the second-order tidal 
approximation, £ cr i t and ji CT n at ISCOs are about 0.0303 and 14.9 irrespective of the value of a (dotted horizontal lines). 
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(c) (d) 

FIG. 5. Density contour curves at Roche limits in the (a) fourth-, (b) third-, and (c) second-order tidal approximations for 
r — 6M, Ro = M, a = 0, and n = 1. (d) density contour curves at a Roche limit in the fourth-order approximation for 
r = 2AM, Ro = M, a = 0.9M , and n = 1. The contour curves are drawn for p/p c = 1CT ' 23 for j = 0, 1, 2, • • • , 15. X and V 
denote x 1 and a; 3 , respectively. 




(a) (b) 

FIG. 6. ("crit and /i C r« as functions of r (a) for a = and (b) for a = 0.9M. For both figures, the results with n = 0.5, 1, and 
1.5 are shown. The solid and dotted curves for each value of n denote the results in the fourth-order tidal approximation with 
Ro = M and in the second-order tidal approximation, respectively. 
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(a) (b) 

FIG. 7. Ccrit and /i cr it as functions of r in the presence of the gravitomagnetic (GM) term (a) for a — and (b) for a = 0.9M. 
For both figures, the results with n = 1 and with Ro/M = 0.5 and 1 are shown. The solid, dashed, long-dashed, and 
dotted-dashed curves denote the results in the fourth-order tidal approximation with the GM term and for Ro = M, in the 
fourth-order tidal approximation with no GM term and for Ro = M, and in the fourth-order tidal approximation with the GM 
term and for Ro = 0.5M, and in the second-order tidal approximation with no GM term, respectively. 




(a) (b) 

FIG. 8. Critical mass M cr it of a black hole for the tidal disruption of a neutron star of mass 1.4M© and radius Ro = 10 km 
and 15 km as a function of r/M (a) for a = and (b) for a — 0.9M. The dashed, solid, and dotted curves denote the results for 
n = 0.5, 1, and 1.5, respectively. For M < M CI i t at a given value of r/M, the neutron star is unstable against tidal disruption. 
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